home *** CD-ROM | disk | FTP | other *** search
/ Whiteline: Alpha / Whiteline Alpha.iso / progtool / modula2 / lpr / longmath.mod < prev    next >
Encoding:
Modula Implementation  |  1994-09-22  |  10.9 KB  |  451 lines
  1. IMPLEMENTATION MODULE LongMathLib0;
  2. FROM Terminal IMPORT WriteString,WriteLn;
  3. FROM SYSTEM   IMPORT LONG;
  4. FROM System   IMPORT FLOATd, TRUNCd;
  5. FROM RealInOut IMPORT Rec, Long;
  6.  
  7. (*      Methods and approximations are taken from 
  8.         HART et al, 
  9.         Computer Approximations, 
  10.         SIAM Series in Applied Mathematics,
  11.         John Wiley and Sons, 
  12.         New York, 1968
  13. *)
  14. TYPE    RealCard = RECORD
  15.                      CASE  : BOOLEAN OF
  16.                        FALSE: x : LONGREAL |
  17.                        TRUE : y : CARDINAL;
  18.                      END;
  19.                    END;
  20.                    
  21. CONST RedMax = 32;
  22.                    
  23. VAR     x,y,z   : LONGREAL;
  24.         NegSign : BOOLEAN;
  25.         shift   : RECORD
  26.                     CASE : BOOLEAN OF
  27.                       TRUE  : I : INTEGER |
  28.                       FALSE : C : CARDINAL
  29.                     END
  30.                   END;
  31.         n       : CARDINAL;
  32.  
  33. (* "CONSTANTS": *)
  34.         Red : ARRAY [1..RedMax] OF LONGREAL;
  35.         Nul,One,Two,Three,Four,Ten,Three70,Half,RTen,TwoPi,PiH,PiQ,Pi2,Pi4,
  36.         TanPi8,Sqrt2,Sqrt21,Ln2,
  37.         EM140,E150,EM150,EM15 : LONGREAL;
  38.         NulI : LONGINT;
  39.  
  40.         SinPar    : RECORD P0,P1,P2,P3,P4,P5,Q0,Q1       : LONGREAL END;
  41.         TanPar    : RECORD P0,P1,P2,P3,P4,   Q0,Q1,Q2    : LONGREAL END;
  42.         LnPar     : RECORD P0,P1,P2,P3,      Q0,Q1,Q2    : LONGREAL END;
  43.         ExpPar    : RECORD P0,P1,            Q0,Q1       : LONGREAL END;
  44.         SqrtPar   : RECORD P0,P1,P2,P3,      Q0,Q1,Q2    : LONGREAL END;
  45.         ArctanPar : RECORD P0,P1,P2,P3,P4,   Q0,Q1,Q2,Q3 : LONGREAL END;
  46.  
  47. PROCEDURE Entier(x : LONGREAL) : LONGINT;
  48.   BEGIN
  49.     NegSign:=x<Nul;
  50.     x:=ABS(x);
  51.     IF x>FLOATd(MAX(LONGINT)) THEN
  52.       WriteLn; 
  53.       WriteString('Error: LONGREAL too large in Entier ');
  54.       WriteLn;
  55.       HALT;
  56.       RETURN 0
  57.     END;
  58.     IF NegSign THEN 
  59.       RETURN -(TRUNCd(x)+LONG(1))
  60.     ELSE 
  61.       RETURN TRUNCd(x) 
  62.     END
  63.   END Entier;
  64.   
  65. PROCEDURE LongReal(x : LONGINT) : LONGREAL;
  66.   BEGIN
  67.     IF x<NulI THEN
  68.       RETURN -FLOATd(ABS(x))
  69.     ELSE
  70.       RETURN FLOATd(x)
  71.     END
  72.   END LongReal;
  73.  
  74. PROCEDURE Ln(A : LONGREAL) : LONGREAL;
  75.   VAR b : RealCard;
  76. (* Hart function 2705 *)
  77.   BEGIN
  78.     IF A<=Nul THEN
  79.       WriteLn;
  80.       WriteString('Error: Ln of value<=0.0 ');
  81.       WriteLn;
  82.       HALT;
  83.       RETURN Nul
  84.     END;
  85.     b.x:=A;
  86.     shift.I:=b.y DIV 10H;
  87.     DEC(shift.I,3FEH);                 (*shift.I by power of 2*)
  88.     b.y:= 3FE0H+(b.y MOD 10H);      (*b now in range 0.5..1.0*)
  89.     z:=b.x;
  90.     WHILE z<Sqrt21 DO 
  91.       DEC(shift.I);
  92.       z:=Two*z;       (*can be improved - use shifts etc*)
  93.     END;
  94.     z:=(z-One)*Rec(z+One);
  95.     y:=z*z;
  96.     WITH LnPar DO
  97.       y:=z*(((P3*y+P2)*y+P1)*y+P0)*Rec(((y+Q2)*y+Q1)*y+Q0)
  98.     END;
  99.         (*      now factor shift.I back in using the fact that 
  100.                 ln(2^n*x)=ln(2^n)+ln(x)=n*ln(2)+ln(x) *)
  101.     RETURN y+LongReal(shift.I)*Ln2
  102.   END Ln;
  103.  
  104. PROCEDURE Exp(A : LONGREAL) : LONGREAL;
  105. (* Hart function 1801*)
  106.   VAR b : RealCard;
  107.       n : INTEGER;
  108.   BEGIN
  109.     IF A<Nul THEN
  110.       b.x:=-A;
  111.       IF b.x>Three70 THEN RETURN Nul END;
  112.       NegSign:=TRUE;
  113.     ELSE
  114.       b.x:=A;
  115.       IF A>Three70 THEN
  116.         WriteLn;
  117.         WriteString('Error: value too large in Exp');
  118.         WriteLn;
  119.         HALT;
  120.         RETURN Nul
  121.       END;
  122.       NegSign:=FALSE;
  123.     END;
  124.     shift.I:=b.y DIV 10H;
  125.     DEC(shift.I,3FBH);             (*shift.I by power of 2*)
  126.     IF shift.I>0 THEN
  127.       b.y:=3FB0H+(b.y MOD 10H)  (*b now in range 0..1/32*)
  128.     ELSE
  129.       shift.I:=0
  130.     END;
  131.     y:=b.x*b.x;
  132.     WITH ExpPar DO
  133.       z:=P1*y+P0;
  134.       y:=One+(Two*b.x*z)*Rec(((y+Q1)*y+Q0)-b.x*z)
  135.     END;
  136.     FOR n:=shift.I TO 1 BY -1 DO
  137.       y:=y*y;
  138.     END;    
  139.     IF NegSign THEN y:=Rec(y) END;
  140.     RETURN y;
  141.   END Exp;
  142.  
  143. PROCEDURE Sqrt(A : LONGREAL) : LONGREAL;
  144.   VAR b   : RealCard;
  145.       exp : CARDINAL;
  146.   BEGIN
  147.     IF A=Nul THEN RETURN A END;
  148.     IF A<Nul THEN
  149.       WriteLn;
  150.       WriteString('Error: Sqrt of value<0.0 ');
  151.       WriteLn;
  152.       HALT;
  153.       RETURN Nul
  154.     END;
  155.     b.x:=A;
  156.     shift.I:=b.y DIV 10H;
  157.     DEC(shift.I,3FEH);          (*shift.I by power of 2*)
  158.     IF ODD(shift.I) THEN 
  159.       exp:=3FD0H;
  160.       INC(shift.I)
  161.     ELSE
  162.       exp:=3FE0H
  163.     END;
  164.     b.y:=exp + b.y MOD 10H;  (*b.x now in [0..1/2]*)
  165.     x:=b.x;
  166.     WITH SqrtPar DO
  167.       y:=((((P3*x)+P2)*x+P1)*x+P0)*Rec(((x+Q2)*x+Q1)*x+Q0)
  168.     END;
  169.         (*this is the first guess at a result*)
  170.     REPEAT
  171.       z:=Half*(x*Rec(y)-y);
  172.       y:=y+z;
  173.     UNTIL ABS(z)<EM15;
  174.     (* now add size back in *)
  175.     IF shift.I<>0 THEN
  176.       b.x:=y;
  177.       INC(b.y,(shift.C DIV 2)*10H);
  178.       RETURN b.x
  179.     ELSE    
  180.       RETURN y;
  181.     END;
  182.   END Sqrt;
  183.   
  184. PROCEDURE Reduce(VAR x : LONGREAL);
  185.   BEGIN
  186.     IF x>Red[RedMax] THEN
  187.       WriteLn;
  188.       WriteString('Error: Too big argument in long trig. function');
  189.       WriteLn;
  190.       HALT;
  191.       x:=Nul;
  192.       RETURN
  193.     END;
  194.     n:=1;
  195.     WHILE x>Red[n] DO INC(n) END;
  196.     IF n>1 THEN
  197.       FOR n:=n-1 TO 1 BY -1 DO
  198.         WHILE x>Red[n] DO x:=x-Red[n] END
  199.       END
  200.     END
  201.   END Reduce;
  202.         
  203. PROCEDURE Sin(A : LONGREAL) : LONGREAL;
  204.   BEGIN
  205.  (* Hart function no 3369 *)
  206.     IF A<Nul THEN
  207.       NegSign:=TRUE;
  208.       x:=-A;
  209.     ELSE
  210.       NegSign:=FALSE;
  211.       x:=A;
  212.     END;
  213.     Reduce(x);
  214.     IF x>Pi THEN
  215.       NegSign:=NOT NegSign;
  216.       x:=x-Pi;
  217.     END;
  218.     IF x>PiH THEN
  219.       x:=Pi-x;
  220.     END;
  221.     x:=Pi2*x;
  222.     y:=x*x;
  223.     WITH SinPar DO
  224.       y:=x*(((((P5*y+P4)*y+P3)*y+P2)*y+P1)*y+P0)*Rec((y+Q1)*y+Q0)
  225.     END;
  226.     IF NegSign THEN 
  227.       RETURN -y
  228.     ELSE 
  229.       RETURN y 
  230.     END;
  231.   END Sin;
  232.  
  233. PROCEDURE Cos(A : LONGREAL) : LONGREAL;
  234.   BEGIN
  235.     RETURN Sin(A+PiH);   (*can do better - eg Hart 3843*) 
  236.   END Cos;
  237.  
  238. PROCEDURE Tan(A : LONGREAL) : LONGREAL;
  239.   PROCEDURE Hart4285(A : LONGREAL) : LONGREAL;
  240.     BEGIN
  241.       x:=Pi4*A;
  242.       y:=x*x;
  243.       WITH TanPar DO
  244.         RETURN x*((((P4*y+P3)*y+P2)*y+P1)*y+P0)*Rec(((y+Q2)*y+Q1)*y+Q0)
  245.       END
  246.     END Hart4285;
  247.   BEGIN
  248.     IF A<Nul THEN
  249.       x:=-A;
  250.       NegSign:=TRUE
  251.     ELSE
  252.       x:=A;
  253.       NegSign:=FALSE
  254.     END;
  255.     Reduce(x);
  256.     IF x>Pi THEN x:=x-Pi END;
  257.     IF x>(PiH) THEN
  258.       x:=Pi-x;
  259.       NegSign:=NOT NegSign;
  260.     END;
  261.     IF x>(PiQ) THEN 
  262.       z:=x;
  263.       x:=(PiH-x);
  264.       IF x<EM140 THEN 
  265.         WriteLn;
  266.         WriteString('Error: Tan too close to Pi/2');
  267.         WriteLn;
  268.         y:=E150;
  269.       ELSE    
  270.         z:=Hart4285(x);
  271.         y:=Rec(z);
  272.       END;
  273.     ELSE 
  274.       y:=Hart4285(x) 
  275.     END;
  276.     IF NegSign THEN 
  277.       RETURN -y
  278.     ELSE 
  279.       RETURN y 
  280.     END;
  281.   END Tan;
  282.  
  283. PROCEDURE Arctan(A : LONGREAL) : LONGREAL;
  284.   VAR NegSign:BOOLEAN;
  285.   PROCEDURE Hart5076(x : LONGREAL) : LONGREAL;
  286.     (*The polynomial approximation to arctan on [0..Pi/8] *)
  287.     BEGIN
  288.       y:=x*x;
  289.       WITH ArctanPar DO
  290.         RETURN x*((((P4*y+P3)*y+P2)*y+P1)*y+P0)*
  291.                Rec((((y+Q3)*y+Q2)*y+Q1)*y+Q0)
  292.       END
  293.     END Hart5076;
  294.   BEGIN
  295.     IF A<Nul THEN
  296.       NegSign:=TRUE;
  297.       A:=-A;
  298.     ELSE
  299.       NegSign:=FALSE;
  300.     END;
  301.     IF A>EM15 THEN 
  302.        (*      Reduce range to be in [0..tan(Pi/8)]*)
  303.       IF A=One THEN 
  304.         A:=PiQ
  305.       ELSIF A>One THEN
  306.         A:=PiH-Arctan(Rec(A))
  307.       ELSIF A>=TanPi8 (*Tan(Pi/8)*) THEN
  308.         A:=PiQ-Hart5076(Two*Rec(One+A)-One)
  309.       ELSE
  310.         A:=Hart5076(A);
  311.       END;
  312.     END;
  313.     IF NegSign THEN 
  314.       RETURN -A
  315.     ELSE 
  316.       RETURN A
  317.     END;
  318.   END Arctan;
  319.  
  320. PROCEDURE Arctan2(Y, X : LONGREAL) : LONGREAL;
  321.   VAR Quadrant  :CARDINAL; 
  322.   BEGIN
  323.     IF (X>=Nul) THEN
  324.       IF Y>=Nul THEN 
  325.         Quadrant:=1
  326.       ELSE
  327.         Y:=-Y;
  328.         Quadrant:=4;
  329.       END
  330.     ELSE
  331.       X:=-X;
  332.       IF Y>=Nul THEN 
  333.         Quadrant:=2
  334.       ELSE
  335.         Y:=-Y;
  336.         Quadrant:=3;
  337.       END;
  338.     END;
  339.     IF X<=EM150*Y THEN 
  340.       x:=PiH
  341.     ELSE x:=Arctan(Y*Rec(X)) END;
  342.     CASE Quadrant OF
  343.       1: RETURN x    |
  344.       2: RETURN Pi-x |
  345.       3: RETURN Pi+x |
  346.       4: RETURN Two*Pi-x
  347.     END;
  348.   END Arctan2;
  349.  
  350. BEGIN (* Initialisation of "Constants" *)
  351.  
  352.   NulI:=0;
  353.   One:=1.0;
  354.   Ten:=(10.0);
  355.   EM15:=1.0E-15;
  356.   EM140:=Long(1000,0,0,0,0,-140);
  357.   E150 :=Long(1000,0,0,0,0, 150);
  358.   EM150:=Long(1000,0,0,0,0,-150);
  359.     
  360.   Nul:=  0.0;
  361.   Two:=  2.0;
  362.   Three:=3.0;
  363.   Four:= 4.0;
  364.   Three70:=370.0;
  365.   Half:=Rec(Two);
  366.   
  367.   Pi:=Long(3141,5926,5358,9793,2300,0);
  368.   TwoPi:=Two*Pi;
  369.   PiH:=Pi*Rec(Two);
  370.   PiQ:=Pi*Rec(Four);
  371.   Pi2:=Two*Rec(Pi);
  372.   Pi4:=Four*Rec(Pi);
  373.   
  374.   Red[1]:=TwoPi;
  375.   FOR n:=2 TO RedMax DO
  376.     Red[n]:=Three*Red[n-1]
  377.   END;
  378.   
  379.   TanPi8:=Long(4142,1356,2373,0950,4880,-1);
  380.   
  381.   E :=Long(2718,2818,2845,9045,2300,0);
  382.  
  383.   Ln2:=Long(6931,4718,0559,9453,0942,-1);
  384.   
  385.   WITH LnPar DO
  386.     P3:= Long(4210,8737,1217,9797,1450,-1);
  387.     P2:=-Long(9637,6909,3368,6865,9324, 0);
  388.     P1:= Long(3095,7292,8215,3765,0062, 1);
  389.     P0:=-Long(2401,3917,9559,2105,0987, 1);
  390.     Q2:=-Long(8911,1090,2793,7831,2337, 0);
  391.     Q1:= Long(1948,0966,0700,8897,3052, 1);
  392.     Q0:=-Long(1200,6958,9779,6052,5472, 1)
  393.   END;
  394.   
  395.   WITH ExpPar DO
  396.     P1:= Long(2000,1114,1589,9645,6894, 1);
  397.     P0:= Long(8400,6685,2536,4832,3941, 2);
  398.     Q1:= Long(1800,1337,0407,3900,2281, 2);
  399.     Q0:= Long(1680,1337,0507,2966,4841, 3)
  400.   END;
  401.   
  402.   WITH SqrtPar DO
  403.     P3:= Long(2975,3039,1000,0000,0000, 0);
  404.     P2:= Long(2027,7246,3000,0000,0000, 0);
  405.     P1:= Long(1095,4240,5000,0000,0000,-1);
  406.     P0:= Long(3162,3500,0000,0000,0000,-4);
  407.     Q2:= Long(3463,9955,6000,0000,0000, 0);
  408.     Q1:= Long(6412,2536,7000,0000,0000,-1);
  409.     Q0:= Long(9408,9090,0000,0000,0000,-3)
  410.   END;
  411.    
  412.   WITH SinPar DO
  413.     P5:=-Long(2276,7796,5988,7676,1976,-3);
  414.     P4:= Long(2573,8460,9831,1601,9960,-1);
  415.     P3:=-Long(1239,4583,0531,8782,6270, 1);
  416.     P2:= Long(2834,8568,2067,7127,2436, 2);
  417.     P1:=-Long(2707,7267,5338,4732,7688, 3);
  418.     P0:= Long(7024,8302,5674,7977,8974, 3);
  419.     Q1:= Long(1153,0387,1236,1456,1194, 2);
  420.     Q0:= Long(4472,1458,3897,1795,8397, 3)
  421.   END;
  422.   
  423.   WITH TanPar DO
  424.     P4:= Long(3386,6386,4267,7172,0961,-5);
  425.     P3:= Long(3422,5543,8724,1003,4353,-2);
  426.     P2:=-Long(1550,6856,5348,3266,3769, 1);
  427.     P1:= Long(1055,9709,0171,4953,1936, 3);
  428.     P0:=-Long(1306,8202,6475,4825,6683, 4);
  429.     Q2:=-Long(1555,0331,6403,1709,9669, 2);
  430.     Q1:= Long(4765,7513,6291,6483,6989, 3);
  431.     Q0:=-Long(1663,8952,3894,7119,0019, 4)
  432.   END;
  433.  
  434.   WITH ArctanPar DO
  435.     P4:= Long(1589,7402,8848,2307,0480,-1);
  436.     P3:= Long(6660,5790,1700,9262,6575, 0);
  437.     P2:= Long(4096,9264,8321,0225,6374, 1);
  438.     P1:= Long(7747,7687,7192,0420,8616, 1);
  439.     P0:= Long(4454,1340,0592,9068,0320, 1);
  440.     Q3:= Long(1550,3977,5514,2198,7525, 1);
  441.     Q2:= Long(6283,5930,5110,3237,6833, 1);
  442.     Q1:= Long(9232,4801,0723,0097,4841, 1);
  443.     Q0:= Long(4454,1340,0592,9068,0445, 1)
  444.   END;
  445.   
  446.   Sqrt2:=Sqrt(Two);
  447.   Sqrt21:=Rec(Sqrt2);
  448.    
  449. END LongMathLib0.
  450.  
  451.